Toward a microscopic description of flow near the jamming 
threshold 
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Abstract -We study the relationship between microscopic structure and viscosity in non- 
Brownian suspensions. We argue that the formation and opening of contacts between particles 
in flow effectively leads to a negative selection of the contacts carrying weak forces. We show 
that an analytically tractable model capturing this negative selection correctly reproduces scaling 
properties of flows near the jamming transition. In particular, we predict that (i) the viscosity 
n diverges with the coordination number z as n ~ (z c — z)~ ( - 3+e ' > ^ 1+e \ (ii) the operator which 
governs flow displays a low- frequency mode that controls the divergence of viscosity, at a frequency 
i^min ~ (z c — z )( 3+0 ' > ^ 2+2e \ and (in) the distribution of forces displays a scale /* that vanishes 
near jamming as /*/(/) ~ (z c — «) 1 ' < - 1+9 ' where characterizes the distribution of contact forces 
P(f) ~ f e at jamming, and where z c is the Maxwell threshold for rigidity. 
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Suspensions are heterogeneous fluids containing solid 
particles, whose viscosity was computed early on by Ein- 
stein [T] and Batchelor [2] in the dilute regime. As the 
packing fraction <f> is increased however, steric hindrance 
becomes dominant and particles move under stress in an 
increasingly coordinated way [SHE]. For non-Brownian 
particles, the viscosity r\ eventually diverges as the sus- 
pension jams into an amorphous solid at some packing 
fraction cf> c . Recently, progress has been made in char- 
acterizing the rheological properties in this limit. It has 
been shown that non-Brownian suspensions [IH9], as well 
as aerial granular flows [TUHH], are characterized by two 
constitutive relations that relate the packing fraction <j> 
and the macroscopic friction /i to the ratio of the local 
shear rate to the pressure. Both constitutive relations 
display singularities as jamming is approached, in par- 
ticular it is observed that r\ ~ (</> c — (f))~ a where a w 2 
[7]. This phcnomcnological description is valid beyond a 
length scale that grows near jamming, and below which 
non-local effects play a role [3ll8lU2] . There is currently no 
accepted microscopic description for this growing length 
scale, nor for the observed constitutive relations. 



Olsson and Teitel [!1[T3] and others [14][l5] have pop- 
ularized a simplified model of non-Brownian suspensions 
where hydrodynamical interactions are neglected. We re- 
fer to the hard-particle limit of this model as the Affinc 



Solvent Model (ASM). The ASM constitutive relations are 
very similar to those found in real suspensions [TB], sup- 
porting that it captures the essential physics near jam- 
ming. In this model it was observed that [16]: (i) there 
exists a scaling relation between the viscosity r\ and the 
coordination z: r\ ~ (z c — z) -2,85 , where z c = 2D [17] and 
D is the spatial dimension, (ii) The dynamics is governed 
by one operator only. The material can thus be character- 
ized by the spectrum of this operator, which contains more 
information than the constitutive relations. In flow, the 
spectrum displays bi-scaling near jamming, with a single 
mode being responsible for the fast divergence of the vis- 
cosity. In this Letter we explain these observations, make a 
new scaling prediction on the distribution of contact forces 
that we confirm empirically, and define and measure two 
new exponents characterizing contact forces both in flow 
and at jamming. 

ASM is fully defined by the following three assumptions: 
(i) hydrodynamic interactions are neglected: the viscous 
drag on a particle is proportional to the difference V n . a . 
between the particle velocity and the imposed velocity of 
the underlying fluid: F = — ^o^n.a. where £0 is a drag co- 
efficient. The flow of the fluid phase is undisturbed by the 
particles and is chosen to be an affine simple shear of strain 
rate 7. (ii) The dynamics is over-damped. (Hi) Particles 
are hard, i.e. cannot overlap, and frictionless. 
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As is also the case for non-Brownian suspensions of 
hard particles [71 IS], within ASM rheological properties 
depend on only one dimensionless parameter, the normal- 
ized pressure p = p p d D ~ 2 /j^Q [IB]. Here p p is the particle 
pressure and d is the particle size. Near jamming the 
particle shear stress a p and pressure p p are proportional 
fx = (Tp/pp—tfi* >0 [TBTTSIITD] . implying that the viscosity 
rj and the renormalized pressure p are proportional too: 
rj = (Tp/'j —> (£,ofi* /d D ~ 2 ) x p. These quantities depend 
only on the geometry of the network formed by particles 
in contact and can be expressed in a compact form |16j : 



n = 
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where Q is the volume of the system. |r) is the vector 
of dimension iV c - the total number of contacts made 
between the N particles - of the distances r a between 
particles in contact. I7) is a vector, also of dimension 
N c , whose components are the variations of the contact 
lengths under an affine simple shear in the (x,y) plane: 
7„ = dr a /dj = (f a - e x ){r a - e y )/r a , where r a is the vector 
between the centers of the particles forming the contact a. 
Af is a symmetric operator of dimension N c x N c , which 
can be written as Af = T^T- T is the operator of dimen- 
sion ND x N c that assigns to any set of N c contact forces 
I/) the associated net unbalanced forces \F) appearing on 
the particles [2D] : 



\F) = T\f) v«, F l = Y,n a J Q 



(3) 



where a,i labels the contacts made by particle i and where 
n a = r a /r a . The non-zero elements of T thus correspond 
to the unit vectors n a . T* is its transpose. 

Previously, we have numerically performed [16] a spec- 
tral analysis oiN = Sw a,2 |/^)(/"l m flow and found that: 
(i) the spectrum of Af displays bi-scaling and consists of 
two structures: one isolated mode of frequency cj m i n , and 
a plateau of modes appearing above a frequency ui* ~ p~ 5 
with S Rj 0.35, as shown in Fig. 2c, d. (it) The divergence 
of the viscosity is governed by the lowest frequency modes 
I /min), which has a finite projection on the shear direction 
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According to Eq. (JlJ this observation implies that 77 ~ p ~ 
l/w min , as shown in Fig. 3c. 

In this Letter we show that these observations, together 
with the dependence of viscosity on coordination, can be 
explained by one assumption only, namely that the config- 
urations with coordination z visited by the dynamics are 
similar to shear-jammed configurations (of friction fi* and 
of coordination z c ) where contacts carrying the smallest 
forces are removed until the coordination is z. The term 
"similarity" is used here to indicate that the rheological 
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Fig. 1: a) Snapshot of a shear-jammed configuration in 2D. 
The thickness of the lines are proportional to the contact 
forces, the components of |/o), see text, b) Distributions of 
the contact forces in shear-jammed configurations; the dashed 
curve is a fit of the 3D distribution to the functional form 
P(x) ~ —1/ log(a;), and the black line is a power-law with ex- 
ponent 9 = 0.18. 



properties of the two models, i.e. the real configurations 
found in flow and the constructed ones, fall in the same 
universality class. The rational for this similarity is that 
during flow, a negative selection of the weak contacts oc- 
cur: indeed only contacts with a small force are fragile, 
i.e. tend to open and disappear as flow progresses. On the 
other hand, new contacts are formed by collision and can 
immediately carry a significant force. Our rule to gener- 
ate configurations is the simplest one that captures such 
a negative selection and can be analyzed analytically. 

To check our assumption, we use an event-driven code 
[21] to simulate flow. We simulate systems under simple 
shear flow with N = 1000 particles in three dimensions, 
using Lees-Edwards periodic boudary conditions [22], at 
the volume fraction <f> = 0.644 < <j) c . Half of the par- 
ticles are small and half are large; we set the diameter 
ratio of small and large particles to be 1.4. At large den- 
sities jamming can occur spontaneously when the coordi- 
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Fig. 2: a) Average spectra D(ui) of the operator Af v.s. frequency to measured for configurations constructed from shear-jammed 
configurations by gradually removing contacts carrying the weakest forces, see text for details. The spectra are measured by 
binning the constructed configurations according to their various pressures p as indicated in the legend. The peak at low 
frequency corresponds to one mode, and has been amplified by some factor for visibility, b) D(uj) v.s. rescaled frequency 
ujp s with 5 — 0.35. Note the collapse of the onset of the plateau of modes, c) D(ui) v.s. frequency ui of the operator Af, 
measured for the configurations generated in simulations of flow [16]. For each pressure, the spectra are averaged over at least 
300 independent flow configurations, d) D(ui) v.s. rescaled frequency up s with 5 — 0.35 [16] . 



nation fluctuates up to z = z c (see [21] for details), gener- 
ating anisotropic configurations with an effective friction 
/i = [i* . We consider shear-jammed configurations that 
jam after a shear strain of at least 200% is imposed on an 
isotropic system. After averaging over 300 shear-jammed 
states we obtain an average friction coefficient (/i*) = 
0.125, with a standard deviation y/ ((/i* — (/i*)) 2 ) = 0.017, 
indicating that finite size fluctuations are small. This re- 
sult is consistent with previous observations [23] showing 
that fi* is a self-averaging quantity, whose standard devi- 
ation decreases as l/\f~N. 

At jamming rj = oo and according to Eq.([T]) there must 
be one normalized mode |/o) such that Af\fo) = 0. Hence- 
forth we use the tilde notation to refer to quantities char- 
acterizing shear-jammed configurations. Thus 

= (/o|Wo) = (Mftflfo) = ||f |/„)|| 2 , (5) 

implying, together with the definition of T, that |/o) is 
the vector of contact forces that maintain force-balance. 
An example of |/o) is shown in Fig.([T]), together with the 
distribution of the contact forces computed over 300 shear- 
jammed configurations. The distribution of low forces is 
of particular importance, and we find that f (//(/)) ~ 
(f/(f)) 6 with 9 « 0.18. This scaling relation holds well 
for two decades for D = 2. Such small exponents are hard 
to distinguish from the case where 8 — with logarith- 
mic correction, although the power law fits data better, 
especially for D = 2, see Fig.([T]). 

We next apply the following procedure to each three- 



dimensional shear-jammed state: contact forces are 
sorted, and the m = N c (z c — z)/z c weakest contacts are 
removed from the contact network, starting from m = 1. 
Physically, removing a contact corresponds to eroding the 
particles at the contact point, so that a finite gap appears. 
Mathematically, one simply does not include these con- 
tacts when computing the operator T defined in Eq.Q. 
For each m we recompute the associated matrix Af = 7 T 
of dimension N c x iV c with N c = N c — m. 

We now argue that these constructed configurations are 
accurate models of configurations actually visited in flow. 
We first compare the spectra of Af in the two cases and 
show that they scale similarly with pressure. Numerical 
diagonalization of Af readily gives the density of states 
D{u>), displayed in Fig[2k for our constructed configura- 
tions and in FigJ^t for configurations from flow simula- 
tions. These quantities are indeed nearly identical, since 
both (i) exhibit a plateau of modes above some frequency 
scale u* ~ p~ s (~ z c — z, see below) as proven by the 
collapse of the plateaus onsets in the right panels of FigEl 
and (m) display a minimum frequency w m i n , which docs 
not scale with pressure as uj* does, as shown by the lack 
of data collapse of the low- frequency peak of D(ui) in the 
right panels of Fig [2] We indeed find that w m i n ~ ^/y/p 
for our constructed configurations shown in Fig[3h, which 
is also true in flow as recalled in Fig|3h- We note that the 
width of the low-frequency peak of D{u) does not grow 
with increasing pressure, which indicates that the fluctua- 
tions in the spectra of Af are indeed well-controlled. This 
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is further supported by the limited spread of the clouds of 
points of Fig|3Ja and c, which represent our entire data- 
set. 

Next, we test whether our scheme correctly predicts the 
divergence of pressure (or viscosity) with coordination ob- 
served in flows. Figj3ja shows the pressure p vs. Sz = z c — z 
for each constructed configuration, and FigJSJi represents 
the same quantity measured in flow. Remarkably, the scal- 
ing law p ~ 5z~ x / s with 5 « 0.35 holds for both ensembles 
of configurations. 

We now perform a scaling analysis of the properties 
of our constructed configurations. We remove a fraction 
8z/z c of the weakest contacts (i.e. the smallest compo- 
nents of |/o)) and consider the operator J\f associated 
with the constructed configuration. We start by deriv- 
ing a upper bound for the minimal eigenvalue u>^ in of 
A/", and an lower bound for the viscosity. M is sym- 
metric hence w^ in < {x\N\x) for any normalized vector 
| as). We consider the vector |/ ), the projection of |/o) on 
the N c — (1 — Sz/z c )N c remaining contacts. Obviously 
lim^oll/oll = ||/ || ee 1, thus: 



< 



(Jo Wo) 
ll/oll 2 




Fig. 3: Minimal frequency cj m i n of the operator M for our con- 
structed configurations a) and for configurations from simula- 
tions of flow c). Pressure p vs. distance to threshold coordina- 
tion Sz = z c — z for constructed configurations b) and for flow 
d) . These data are well-captured by the relations u m i n ~ 1/ ,/p 
and p ~ Sz^ 1 / 6 with 8 ta 0.35 for both ensembles of configura- 
tions. 



:(/ |AA|/ )^||T|/o)|| 2 = ||F || 2 , (6) 



where ||Fb|| is the norm of the unbalanced force field as- 
sociated with the contact force |/o). Removed contacts 
from a jammed configuration leads to unbalanced forces 
on particles which lost one or more contacts. In partic- 
ular, when a contact a with a force f a is removed, each 
of the two adjacent particles carries an unbalanced total 
force of amplitude f a . In the limit Sz << 1 most parti- 
cles, whose total force is unbalanced, have only lost one 
contact. Thus: 



\Fn 



2^jt =2JV C f 
Jo 



P(f)fdf, 



(7) 



where the sum is taken over all the removed contacts a r . 
Here P(f) is the distribution of the components of the 
vector |/o) and /* is defined as J Q P(f)df = 8z/z c . As 
FigQ] indicates, we observe that at low forces, P(f/(f)) ~ 
(fVN c ) e , where the factor ~ (/) stems from the 

normalization of the vector |/o). We thus obtain two of 
our main predictions: 



if) 



< 



8z 1 + e . 



3 + 9 

Szi+z . 



(8) 
(9) 



Eq.© predicts the emergence of a characteristic force 
scale in flows, that vanishes in relative terms as jamming 
is approached. In order to test this prediction we measure 
in flow the distributions of contact forces for various pres- 
sures, see FigHk, We observe an erosion in the distribution 
P(f/p) of relative contact forces f/p, with a characteristic 
relative force that indeed decays near jamming. To probe 



the scaling of this characteristic force, we seek a rescal- 
ing of the contact force ///* by some /* that collapses 
the low-force tail of the distribution. The best collapse is 
found for /* ~ p5z ~ p x ~ s . As indicated in Table 1, this 
finding corresponds to our prediction for 9 — 0, and is still 
very close to our prediction using 8 = 0.18. We note the 
difficulty in extracting the force scale /*, as the crossover 
of the distributions of forces towards the eroded regime 
is rather weak. Interestingly, we find that the rescaled 
distributions P(f/f*) scale as (/ '/ f*) x with an exponent 
X 0.38. 

We now test the inequality ((9]). Assuming that this 
upper bound is saturated, as excepted if the present vari- 
ational argument captures the essential physics, leads to 
a prediction for the scaling relation between w m i n and 5z. 
As shown in Table 1, this prediction is in very good agree- 
ment with our observations. 

As noted above, D(uj) also displays a frequency scale 
uj* ~ Sz >> cJ m in above which a plateau of modes appears. 
For completeness, we comment on D(u>) in the frequency 
range [w m j n ,o;*]. Although it cannot be observed in our 
numerics due to the limited size of our systems, normal 
modes must be present in this interval. Indeed, a local 
operator like J\f cannot have a single lowest eigenvalue 
w^jjj separated from the rest of the spectrum. This can 
be seen as follows: if |/ m i n ) is the eigenvector correspond- 
ing to the eigenvalue w 2 lin , consider the family of eigen- 
vectors build by modulating |/ m i n ) by plane waves: 
fq,a = exp(i(f • a)/ m in,Q, where a is the position of contact 
a. The set of vectors \f^} are approximatively orthonor- 
mal (fq>\fq) ~ fi^qf Furthermore, for small wave vectors 
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Fig. 4: a) Distributions P(f/p) of contact forces / rescaled 
by the pressure p measured in flow and shear-jammed con- 
figurations. Measurements in flow are done at various p, as 
indicated in the legend. P(f/p) erodes at low f/p as the pres- 
sure diminishes away from the jamming threshold, b) Distri- 
butions of contact forces vs. rescaled force f/(pSz) ~ / '/p 1 ^ 1 '. 
This rescaling leads to the collapse of the distribution at small 
force, indicating the presence of a scaling law. 

the force balance that was nearly satisfied in |/ m i n ) is only 
weakly perturbed, and one finds after a simple calculation 
that ujl = (f$\N\fq) = w 2 lin + Bq 2 , where B is a con- 



stant of order one. Using that the density of wave- vectors 
grows as D(q) ~ q ^ 1 , the density of states of the frequen- 
cies Wg- must then grow as D(uj) oc uj(uj 2 — ijj 2 ^ n )( D ~ 2 ^ 2 . 
Although the u'i are not exact eigenvalues of A/", we ex- 
pect our estimate for D(u) to be qualitatively correct for 

We now derive the divergence of the viscosity 77 (or 
equivalently p since rj ~ p) with coordination. Using the 
convexity of the inverse function and Eq.® leads to: 

(/olATU) > 77TT7|7T ~ Sz~^ (10) 
\jo\JV I Jo) 

The finite friction coefficient /1* in shear-jammed configu- 
rations implies that contact forces have a finite projection 
onto shear: 



Hm^ (7l/o)/||7ll = (7l/o)/||7ll=C>0 
Together with Eas- flJID]) . Eq. (fTTj) leads to: 
V (7l-AT- 1 | 7 )^(7l/o x2 



in 



CI 







-(/olAT-Vo) >C 2 ^-Sz 



fi 



(11) 



3 + 9 



(12) 

where we have neglected the contribution to the viscosity 
stemming from the components of I7) orthogonal to |/o), 
which we expect to be small. From the definition of I7) it is 
clear that ||7|| 2 ~ d 2 N c , whereas fl ~ N c d D . Eq.$T2§ thus 
yields the following lower bound for the rescaled viscosity 
r h , = V /(d D -%): 



fir 



_ 3+9 
> ASz T+TT 



(13) 



where A is a constant. Using that p ~ r\ r we can di- 
rectly compare this prediction with the observations of 
Fig-©- Table 1 shows that the observed scaling law is 
well-captured by the saturation of Eq. (fT3|) . 
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0.84 
2.69 
-2.69 


1 

3 
-3 


1 

2.85 
-2.85 



Table 1: Predicted scaling exponents and exponents measured 
numerically in flow. Comparisons are made using the best fit- 
ted exponent 8 = 0.18 of the force distribution P{f) of jammed 
states shown in FigfJ] and 9 — assuming that the decay in 
P(f) is logarithmic. 



Conclusion. — In flow, the contact network con- 
stantly evolves, or rewires, via the formation and open- 
ing of contacts. We have shown that a key aspect of this 
process is the negative selection of contacts that weakly 
affect flow. Taking this effect into account enables one 
to derive three scaling relations between four exponents. 
These relations connect the divergence of the viscosity and 
spectral properties of dense flows to two microscopic quan- 
tities: the coordination z and the exponent 9, which char- 
acterizes the density of weak contact forces in jammed 
configurations. Thus to obtain a complete description of 
the rheology, future works should compute the value of 9 
- see [24] for recent results in this direction- and the rela- 
tion between coordination and packing fraction z{<p). In 
packings of soft repulsive particles with <fi > <fi c , the co- 
ordination is the minimal one that guarantees mechanical 
stability a condition that determines z(qA. The 

concept of stability is not applicable for fluids however, 
and computing z(<fi) in flow will presumably require a de- 
tailed description of the rewiring dynamics. 
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